home *** CD-ROM | disk | FTP | other *** search
/ Almathera Ten Pack 3: CDPD 3 / Almathera Ten on Ten - Disc 3: CDPD3.iso / fish / 701-725 / 709 / planets / moon_pos.c < prev    next >
C/C++ Source or Header  |  1995-03-18  |  13KB  |  246 lines
  1. /*******************************************************************************
  2. ** MOON-POS.C - Determine the position of the moon.                           **
  3. **                                                                            **
  4. ** Author: Jim Cobb                                                           **
  5. **         jcobb@utah.edu.cs                                                  **
  6. **         Computer Science Dept.                                             **
  7. **         University of Utah                                                 **
  8. ** Date:   Fri Mar 25 1988                                                    **
  9. **                                                                            **
  10. ** Copyright (c) 1988, Jim Cobb, All rights reserved.                         **
  11. **                                                                            **
  12. ** (This copyright prevents you from selling the program - the author grants  **
  13. ** anyone the right to copy and install the program on any machine it will    **
  14. ** run on)                                                                    **
  15. **                                                                            **
  16. ** Reference: Jean Meeus, "Astronomical Formulae for Calculators," chapter 30 **
  17. ** Meeus states that this procedure has an accuracy of 10" in the longitude   **
  18. ** of the moon, 3" in its latitude and 0.2" in its parallax.                  **
  19. *******************************************************************************/
  20.  
  21. #include <stdio.h>
  22. #include <math.h>
  23.  
  24. #ifndef PI
  25. #  define PI 3.14159265358979323846 /* Should be in <math.h> */
  26. #endif /* PI */
  27. #ifndef RAD
  28. #  define RAD 0.01745329251994329651 /* (PI / 180 deg) */
  29. #endif /* RAD */
  30.  
  31. /*******************************************************************************
  32. ** Trigonometric functions in degrees.                                        **
  33. *******************************************************************************/
  34. #define sind(x) sin(RAD*(x))
  35. #define cosd(x) cos(RAD*(x))
  36.  
  37. #define MAGMOON -9.9
  38.  
  39. /*******************************************************************************
  40. ** "moon_pos" takes a time argument, t, the number of Julian centuries from   **
  41. ** 1900 January 0.5 ET, and the obliquity of the ecliptic, epli.              **
  42. ** It adds information about the moon's position to the logfile.              **
  43. *******************************************************************************/
  44. void moon_pos(t,epli)
  45. double t,epli;
  46.    {
  47.    double e,e_2,
  48.           m,             /* Sun's anomaly.                              */
  49.           m_prime,       /* Moon's anomaly.                             */
  50.           d,             /* Moon's elongation.                          */
  51.           f,             /* Distance of Moon from its ascending node.   */
  52.           omega,         /* Longitude of Moon's ascending node.         */
  53.           venus_term,    /* Great Venus term.                           */
  54.           lambda,        /* Ecliptical longitude of the Moon's center.  */
  55.           beta;          /* Ecliptical latitude of the Moon's center.   */
  56.  
  57.    double phase; /* Phase of the moon */
  58.  
  59. #ifdef EUTSCH
  60.    static char namestring[] = "Mond, 180.";
  61. #else
  62.    static char namestring[] = "Moon, 180.";
  63. #endif
  64.  
  65.    double range(),poly();
  66.    int    to_int();
  67.    void   geo_trans();
  68.  
  69. /*******************************************************************************
  70. ** Compute the mean values for the terms.                                     **
  71. *******************************************************************************/
  72.    lambda  = poly(270.434164,481267.8831,-0.001133,-0.0000019,t);
  73.    m       = poly(358.475833, 35999.0498,-0.000150, 0.0000033,t);
  74.    m_prime = poly(296.104608,477198.8491, 0.009192, 0.0000144,t);
  75.    d       = poly(350.737486,445267.1142,-0.001436,-0.0000019,t);
  76.    f       = poly( 11.250889,483202.0251,-0.003211, 0.0000003,t);
  77.    omega   = poly(259.183275, -1934.1420,-0.002078,-0.0000022,t);
  78.  
  79. /*******************************************************************************
  80. ** Additive terms.                                                            **
  81. *******************************************************************************/
  82.    e_2      = sind( 51.2 + 20.2 * t );
  83.    lambda  +=  0.000233 * e_2;
  84.    m       += -0.001778 * e_2;
  85.    m_prime +=  0.000817 * e_2;
  86.    d       +=  0.002011 * e_2;
  87.  
  88.    venus_term  = 0.003964 * sind(poly(346.560,132.870,-0.0091731,0.,t));
  89.    lambda     += venus_term;
  90.    m_prime    += venus_term;
  91.    d          += venus_term;
  92.    f          += venus_term;
  93.  
  94.    e_2      = sind( omega );
  95.    lambda  +=  0.001964 * e_2;
  96.    m_prime +=  0.002541 * e_2;
  97.    d       +=  0.001964 * e_2;
  98.    f       += -0.024691 * e_2;
  99.  
  100.    f += -0.004328 * sind( omega + 275.05 - 2.30 * t );
  101.  
  102.    e   = poly(1.,-0.002495,-0.00000752,0.,t);
  103.    e_2 = e*e;
  104.  
  105. /*******************************************************************************
  106. ** Bring these angles within 0 to 360 degrees.                                **
  107. *******************************************************************************/
  108.    m       = range( m );
  109.    m_prime = range( m_prime );
  110.    d       = range( d );
  111.    f       = range( f );
  112.    omega   = range( omega );
  113.  
  114. /*******************************************************************************
  115. ** Calculate lambda, the ecliptical longitude of the Moon's center.           **
  116. *******************************************************************************/
  117.    lambda +=  6.288750 * sind(            m_prime               )
  118.       +       1.274018 * sind(  2.*d -    m_prime               )
  119.       +       0.658309 * sind(  2.*d                            )
  120.       +       0.213616 * sind(         2.*m_prime               )
  121.       - e   * 0.185596 * sind(                         m        )
  122.       -       0.114336 * sind(                             2.*f )
  123.       +       0.058793 * sind(  2.*d - 2.*m_prime               )
  124.       + e   * 0.057212 * sind(  2.*d -    m_prime -    m        )
  125.       +       0.053320 * sind(  2.*d +    m_prime               )
  126.       + e   * 0.045874 * sind(  2.*d              -    m        )
  127.       + e   * 0.041024 * sind(            m_prime -    m        )
  128.       -       0.034718 * sind(     d                            )
  129.       - e   * 0.030465 * sind(            m_prime +    m        )
  130.       +       0.015326 * sind(  2.*d                     - 2.*f )
  131.       -       0.012528 * sind(            m_prime        + 2.*f )
  132.       -       0.010980 * sind(       -    m_prime        + 2.*f )
  133.       +       0.010674 * sind(  4.*d -    m_prime               )
  134.       +       0.010034 * sind(         3.*m_prime               )
  135.       +       0.008548 * sind(  4.*d - 2.*m_prime               )
  136.       - e   * 0.007910 * sind(  2.*d -    m_prime +    m        )
  137.       - e   * 0.006783 * sind(  2.*d              +    m        )
  138.       +       0.005162 * sind(  -  d +    m_prime               )
  139.       + e   * 0.005000 * sind(     d              +    m        )
  140.       + e   * 0.004049 * sind(  2.*d +    m_prime -    m        )
  141.       +       0.003996 * sind(  2.*d + 2.*m_prime               )
  142.       +       0.003862 * sind(  4.*d                            )
  143.       +       0.003665 * sind(  2.*d - 3.*m_prime               )
  144.       + e   * 0.002695 * sind(         2.*m_prime -    m        )
  145.       +       0.002602 * sind( -2.*d +    m_prime        - 2.*f ) 
  146.       + e   * 0.002396 * sind(  2.*d - 2.*m_prime -    m        )
  147.       -       0.002349 * sind(     d +    m_prime               )
  148.       + e_2 * 0.002249 * sind(  2.*d              - 2.*m        )
  149.       - e   * 0.002125 * sind(         2.*m_prime +    m        )
  150.       - e_2 * 0.002079 * sind(                      2.*m        )
  151.       + e_2 * 0.002059 * sind(  2.*d -    m_prime - 2.*m        )
  152.       -       0.001773 * sind(  2.*d +    m_prime        - 2.*f )
  153.       -       0.001595 * sind(  2.*d                     + 2.*f )
  154.       +  e  * 0.001220 * sind(  4.*d -    m_prime -    m        )
  155.       -       0.001110 * sind(         2.*m_prime        + 2.*f )
  156.       +       0.000892 * sind( -3.*d +    m_prime               )
  157.       - e   * 0.000811 * sind(  2.*d +    m_prime +    m        )
  158.       + e   * 0.000761 * sind(  4.*d - 2.*m_prime -    m        )
  159.       + e_2 * 0.000717 * sind(            m_prime - 2.*m        )
  160.       + e_2 * 0.000704 * sind( -2.*d +    m_prime - 2.*m        )
  161.       + e   * 0.000693 * sind(  2.*d - 2.*m_prime +    m        )
  162.       + e   * 0.000598 * sind(  2.*d              -    m - 2.*f )
  163.       +       0.000550 * sind(  4.*d +    m_prime               )
  164.       +       0.000538 * sind(         4.*m_prime               )
  165.       + e   * 0.000521 * sind(  4.*d              -    m        )
  166.       +       0.000486 * sind( -   d + 2.*m_prime               );
  167.  
  168.    lambda = range( lambda );
  169.  
  170. /*******************************************************************************
  171. ** Calculate beta, the ecliptical latitude of the Moon's center.              **
  172. *******************************************************************************/
  173.    beta =     5.128189 * sind(                                f )
  174.       +       0.280606 * sind(            m_prime        +    f )
  175.       +       0.277693 * sind(            m_prime        -    f )
  176.       +       0.173238 * sind(  2.*d                     -    f )
  177.       +       0.055413 * sind(  2.*d -    m_prime        +    f )
  178.       +       0.046272 * sind(  2.*d -    m_prime        -    f )
  179.       +       0.032573 * sind(  2.*d                     +    f )
  180.       +       0.017198 * sind(         2.*m_prime        +    f )
  181.       +       0.009267 * sind(  2.*d +    m_prime        -    f )
  182.       +       0.008823 * sind(         2.*m_prime        -    f )
  183.       + e   * 0.008247 * sind(  2.*d              -    m -    f )
  184.       +       0.004323 * sind(  2.*d - 2.*m_prime        -    f )
  185.       +       0.004200 * sind(  2.*d +    m_prime        +    f )
  186.       + e   * 0.003372 * sind( -2.*d              -    m +    f )
  187.       + e   * 0.002472 * sind(  2.*d -    m_prime -    m +    f )
  188.       + e   * 0.002222 * sind(  2.*d              -    m +    f )
  189.       + e   * 0.002072 * sind(  2.*d -    m_prime -    m -    f )
  190.       + e   * 0.001877 * sind(            m_prime -    m +    f )
  191.       +       0.001828 * sind(  4.*d -    m_prime        -    f )
  192.       - e   * 0.001803 * sind(                         m +    f )
  193.       -       0.001750 * sind(                             3.*f )
  194.       + e   * 0.001570 * sind(            m_prime -    m -    f )
  195.       -       0.001487 * sind(     d                     +    f )
  196.       - e   * 0.001481 * sind(            m_prime +    m +    f )
  197.       + e   * 0.001417 * sind(       -    m_prime -    m +    f )
  198.       + e   * 0.001350 * sind(                    -    m +    f )
  199.       +       0.001330 * sind( -   d                     +    f )
  200.       +       0.001106 * sind(         3.*m_prime        +    f )
  201.       +       0.001020 * sind(  4.*d                     -    f )
  202.       +       0.000833 * sind(  4.*d -    m_prime        +    f )
  203.       +       0.000781 * sind(            m_prime        - 3.*f )
  204.       +       0.000670 * sind(  4.*d - 2.*m_prime        +    f )
  205.       +       0.000606 * sind(  2.*d                     - 3.*f )
  206.       +       0.000597 * sind(  2.*d + 2.*m_prime        -    f )
  207.       + e   * 0.000492 * sind(  2.*d +    m_prime -    m -    f )
  208.       +       0.000450 * sind( -2.*d + 2.*m_prime        -    f )
  209.       +       0.000439 * sind(         3.*m_prime        -    f )
  210.       +       0.000423 * sind(  2.*d + 2.*m_prime        +    f )
  211.       +       0.000422 * sind(  2.*d - 3.*m_prime        -    f )
  212.       - e   * 0.000367 * sind(  2.*d -    m_prime +    m +    f )
  213.       - e   * 0.000353 * sind(  2.*d              +    m +    f )
  214.       +       0.000331 * sind(  4.*d                     +    f )
  215.       + e   * 0.000317 * sind(  2.*d +    m_prime -    m +    f )
  216.       + e_2 * 0.000306 * sind(  2.*d              - 2.*m -    f )
  217.       -       0.000283 * sind(            m_prime        + 3.*f );
  218.  
  219.    beta *= ( 1.
  220.       - 0.0004664 * cosd( omega )
  221.       - 0.0000754 * cosd( omega + 275.05 - 2.30 * t ));
  222.  
  223. /*******************************************************************************
  224. ** Roughly calculate the phase of the moon.                                   **
  225. ** Added by ccount                                                            **
  226. ** Sun longitude is about 36000 * T0 + 279                                    **
  227. ** Phase is 0 for full moon                                                   **
  228. *******************************************************************************/
  229.     e_2   = 36000. * t + 279 - lambda;
  230.     phase = 180. - (e_2 - ((double)(to_int(e_2/360.))*360.));
  231.  
  232.     if(phase < 0.)
  233.        phase += 360.;
  234.  
  235. #ifdef EUTSCH
  236.     sprintf(namestring,"Mond, %3.0f",phase);
  237. #else
  238.     sprintf(namestring,"Moon, %3.0f",phase);
  239. #endif
  240.  
  241. /*******************************************************************************
  242. ** Transform to right ascension and declination, and output the result.       **
  243. *******************************************************************************/
  244.    geo_trans(lambda,beta,epli,0.,MAGMOON,"PL",namestring );
  245.    }
  246.